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We study Bose-Einstein condensation (BEC) in three-dimensional two-component bosonic gases, 
characterizing the universal behaviors of the critical modes arising at the BEC transitions. For 
this purpose, we use field-theoretical (FT) renormalization-group (RG) methods and perform mean- 
field and numerical calculations. The FT RG analysis is based on the Landau-Ginzburg-Wilson 
<E>"^ theory with two complex scalar fields which has the same symmetry as the bosonic system. In 
particular, for identical bosons with exchange Z 2 symmetry, coupled by effective density-density 
interactions, the global symmetry isZ2,e0U(l)(8)U(l). At the BEC transition it may break 
into Z2,e 0 Z2 0 Z2 when both components condense simultaneously, or to U(l) 0 Z2 when only 
one component condenses. This implies different universality classes for the corresponding critical 
behaviors. Numerical simulations of the two-component Bose-Hubbard model in the hard-core limit 
support the RG prediction: when both components condense simultaneously, the critical behavior is 
controlled by a decoupled XY fixed point, with unusual slowly-decaying scaling corrections arising 
from the on-site inter-species interaction. 

PACS numbers: 67.25.dj,67.85.Hj,05.70.Jk,05.10.Cc 


I. INTRODUCTION 

Experiments with cold atomsid— have provided the 
opportunity to investigate Bose-Einstein condensation 
(BEC) in dilute interacting atomic gases. In the BEC 
a macroscopic number of bosonic atoms, the so-called 
condensate, occupy the lowest-energy quantum state at 
a finite temperature. The phase of the condensate wave 
function provides the order parameter at the transition. 
BEC transitions are generically expected to belong to the 
three-dimensional (3D) XY universality class, which is 
characterized by the spontaneous breaking of an Abelian 
U(l) symmetry. The same universal critical behavior is 
observed in the superfiuid transition in in transi¬ 

tions characterized by density or spin waves (as it occurs 
in some liquid crystals), in magnetic systems with easy- 
plane anisotropy, etci^ The XY behavior at the 3D BEC 
transition has been supported by experimental measure¬ 
ments of the diverging correlation length in a cold-atom 
bosonic gasi^ Cold-atom experiments have been extended 
to mixtures of homonuclear and heteronuclear bosonic 
gaseS)^^— which also show BEC phenomena. Several the¬ 
oretical studies have discussed various aspects of the be¬ 
havior of mixtures of bosonics gases, see, e.g., Refs. [ 13 - 
such as low-dimensional behaviors, magnetic-like be¬ 
haviors at the n = 1 Mott phases, etc. However, some 
issues call for further investigations, such as the critical 
behaviors at the finite-temperature normal-to-superfiuid 
transitions which arise from different BEC patterns. 

In this paper we study BEC in mixture of 3D bosonic 
gases, focussing on the critical behaviors at the finite- 
temperature transitions arising from BEC. In particular, 
we consider a system of two identical boson gases with 
density-density interactions. Equivalently, we may inter¬ 
pret this system as made up by a single two-component 
boson gas. An example is provided by the lattice two- 


component Bose-Hubbard (2BH) model 

H = -t^^{blj,y+h.c) - ( 1 ) 

s (xy) sx 

4” 2^ ^ ^ 1) “h C/ ^ ^ f^lx‘^2x') 

sx X 

where (xy) indicates the nearest-neighbor sites of a cu¬ 
bic lattice, the subscript s labels the two species, and 
nsx = is the density operator. The Hamiltonian 

is symmetric under U(l) transformations acting indepen¬ 
dently on the two species and under the Z 2 transforma¬ 
tion exchanging the two bosons. The two-component bo¬ 
son gas shows a quite complex phase diagram in the space 
of the model parameters, i.e., the temperature T, the 
chemical potential fj,, and the on-site couplings U and V. 
In the following we set t = 1 for the hopping parameter 
without loss of generality. 

We investigate the critical behavior of systems like 
the 2BH model by field-theoretical (FT) renormalization- 
group (RG) methods, mean-field and numerical ap¬ 
proaches. We show that transitions in these two- 
component systems may be associated with different 
spontaneous breakings of the global symmetry 

Z2.e0U(l)0U(l). (2) 

This symmetry may break to 112 ,e 0 Z 2 0 Z 2 when both 
components condense simultaneously, or to U(l) 0 Z 2 
when only one component condenses, with two different 
universality classes for the corresponding critical behav¬ 
iors. 

When both components condense simultaneously, the 
RG analysis shows that the critical behavior is controlled 
by a decoupled 3D XY fixed point (FP). Thus, the tran¬ 
sition belongs to the 3D XY universality class associ¬ 
ated with the symmetry breaking U(l)—>-22. However, 
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the irrelevant density-density interaction between the two 
components gives rise to scaling corrections that decay 
very slowly, as ^-0 022 ^ where ^ is the diverging length 
scale at the transition. Such scaling corrections are not 
present in standard transitions belonging to the XY uni¬ 
versality class, such as at the BEC transition of a single 
bosonic species4^“— In that case, scaling corrections de¬ 
crease significantly faster, as ^ “0.78^ instead, only one 
component condenses, the RG analysis predicts a differ¬ 
ent critical behavior, which belongs to the same univer¬ 
sality class as that of the continuous transitions in chiral 
models with 0(2)(8>0(2) symmetry4S We also present a 
RG analysis of the behavior of mixtures of nonidentical 
bosons with density-density interactions. In this case the 
phase diagram is characterized by the presence of bicrit- 
ical or tetracritical points where various transition lines 
meet. 

The paper is organized as follows. In Sec.HIlwe present 
the RG analysis of the Landau-Ginzburg-Wilson (LGW) 
theory, which is expected to describe the critical tran¬ 
sitions in two-component bosonic systems with density- 
density interactions. In Sec. Hill we discuss the phase 
diagram of the 2BH model (H} in the mean-field approx¬ 
imation. Sec. HYl is devoted to a numerical study of the 
2BH model in the hard-core, R —)> 00 , limit. At the tran¬ 
sition both components condense. We show that the re¬ 
sults can be explained by a 3D XY critical behavior with 
slowly-decaying scaling corrections, as predicted by the 
RG analysis. Finally, in Sec. we draw our conclusions. 


II. FIELD-THEORETICAL 
RENORMALIZATION-GROUP ANALYSIS 

We wish now to classify the finite-temperature transi¬ 
tions in the phase diagram of systems consisting of two 
identical boson species with density-density interactions. 
For this purpose, we study the RG flow of the effective 
LGW theory associated with the critical modesi^i^^— 
Within the FT RG approach, one first identifies the or¬ 
der parameter. Then, one considers the <I>^ Hamiltonian 
with the most general fourth-order potential in the or¬ 
der parameter that has the same symmetry properties 
as the original system. The possible critical behaviors 
are determined by the stable FPs of the RG flow. Each 
of them corresponds to a different universality class, as¬ 
sociated with the symmetry breaking that occurs in the 
parameter region in which the FP is located. The stable 
FPs determine the universal scaling properties, such as 
the critical exponents, the scaling functions, etc. Note 
that only systems which are in the attraction domain 
of the stable FPs undergo continuous transitions. Sys¬ 
tems corresponding to LGW theories with parameters 
that are outside the FP attraction domains, or that be¬ 
long to the instability region, are predicted to undergo 
first-order transitions. 

At a BEC transition, the condensate behaves like the 
magnetization in magnetic systems, i.e., {bax) ~ (Tc ~ 


T)P as Tc is approached from below. Critical modes 
develop a diverging length scale ^ ~ |T — while 

the two-point function at the critical point decays al¬ 
gebraically as^i^ G{x) ^ The exponents /3, 

V, and rj are universal—they only depend on the uni¬ 
versality class—and are related by the scaling relation 
13 = iy{l + T])/2. 


A. LGW theory for two-component boson gases 

In the case of a mixture of bosonic gases, we associate a 
complex field ips{x), s = 1,2, with each bosonic species. 
Since we consider finite-temperature transitions of 3D 
quantum systems, we must consider a three-dimensional 
LGW model. As mentioned in the introduction, the rel¬ 
evant symmetry of the systems we consider, such as the 
2BH model H]), is 1^2,e ® U(l) (Si U(l). The Hamiltonian 
is therefore 

"Hlgw = I |(ps(a:)P (3) 

S,/J. S 

S 

where the potential is the most general one under sym¬ 
metry The Hamiltonian is bounded from below for 
g > 0 and g + u > 0. The quartic couplings g and u are 
related to the intra-species V and inter-species U on-site 
couplings of the 2BH model ([T]). In particular, u must 
vanish when U vanishes, leaving two decoupled LGW 
theories, one for each bosonic gas. 

General information on the phase diagram of model ([3]) 
can be inferred by a straightforward mean-field analysis, 
e.g., by determining the minima of the potential 

V{(p) = + g'^\^s\^+ 2u\ipi\^\ip2\^. (4) 

s s 

For r > 0 the potential is minimized by (ps = 0, while for 
r < 0 the minimum depends on the sign oi w = g — u. If 
w > 0, the minimum occurs when both field components 
condense, i.e., for {(pi) = ((^ 2 ) ^ 0. This implies the 
symmetry-breaking pattern 

^2,6 G u(l) 0 u(l) —>■ 1j2,e 0 Z 2 0 ^2, (5) 

i.e., each U(l) group breaks into Z 2 . Instead for w < 0, 
we have {(pi) 0 and {(P 2 ) = 0, or viceversa. Thus, the 

exchange symmetry and only one of the two U(l) groups 
are broken, so that 

Z2.e®U(l)(g)U(l)^U(l)(g)Z2. (6) 

On the boundary line w = 0, the LGW theory ([3]) is 
equivalent to the 0(4) vector model. 
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FIG. 1: Sketch of the RG flow of the MN model (O for 
M = N = 2. The relevant stable FPs are the decoupled XY 
FP (dXY), corresponding to FP (iii) in the list reported in 
the text, and FP M, corresponding to FP (iv). 


B. RG flow and critical behaviors 


FIG. 2: Sketch of the RG flow of the LGW theory for a mix¬ 
ture of two identical bosonic gases. The relevant stable FPs 
are the decoupled XY FP (dXY), which controls the RG flux 
for ^ — u > 0 (in this case the exchange symmetry is con¬ 
served), and a second FP (Asy), relevant for g — u < 0 (the 
exchange symmetry is broken). 


The LGW theory ([3]) is a particular case of the so- 
called MN model^i^ 


as follows: 


^11 



where (j)ai is an M x N matrix, i.e., a = 1,... ,N and 
i = 1,..., M. Indeed, Hamiltonian ([3]) reduces to ([7]) for 
M = iV = 2, if we set (f)ii = Re(/3i, (j) 2 i = Im(/?i, and 


(/'12 

4>21 

4>22 


- Imv32), 

- Rev?2), 

-^(Im(/3i -I- Rev?2), 

^(Re:/Ji -b lmip2), 
v2 


( 10 ) 


g = vi+V 2 , u = vi, w = g — u = V 2 - (8) 


The RG flow of the MN models has been studied by 
various FT methodsi^ i^°i^^ A sketch of the RG flow in the 
case M = N = 2 is shown Fig. [TJ There are several FPs 
in the plane of the renormalized quartic couplings vi and 
V 2 '^ (i) the trivial Gaussian FP for vi = V 2 = 0 which 
is unstable against both quartic perturbations present in 
Hamiltonian ©; (ii) the 0(4)-symmetric FP for r;2 = 0 
and > 0 which is unstable with respect to the quartic 
term proportional to V 2 in Eq. ©; (iii) a stable decoupled 
XY FP with ui = 0 and V 2 > 0, with attraction domain 
in the region V 2 > 0; (iv) a stable FP for z;2 < 0 with 
attraction domain in the region ^2 < 0. 

It is also possible to map the LGW theory ([3]) onto the 
chiral LGW theory with 0(iV)G>0(M) symmetry defined 

p ..49.54.55 


/ 


d^x 


+ {uo - Vo) + 

ai aij 


VQ 




abij 


'I 


(9) 


where (j)ai is an MxN matrix. The two models are equiv¬ 
alent ioi M = N = 2, if we identify helds and couplings 


and g = uq — vol2., u = uq + uo/2. Note that the cou¬ 
plings of the chiral and of the MN model are related by 
Vo = —V 2 and ^0=1^1+ ^2/2, so that the FP (iv) men¬ 
tioned above corresponds to the chiral FP that has been 
extensively discussed in Refs. 

According to the correspondence ([8]), the equivalent 
model m presents two stable FPs with attraction do¬ 
mains separated by the line w = 5 — u = 0, along which 
the unstable 0(4)-symmetric FP is located. The corre¬ 
sponding RG flows are sketched in Fig. [2j 

If ui > 0 the finite-temperature transition is character¬ 
ized by the simultaneous condensation of both species. 
The stable FP which determines the critical behavior 
is located along the line u = 0, thus representing two 
decoupled U(l)-symmetric models. This implies that it 
belongs to the 3D XY universality class, whose critical 
exponents are^ Vxy = 0.6717(1) and tjxy = 0.0381(2). 
However, scaling corrections decay much slower than in 
the U(l)-symmetric theory with a single complex field, 
where the leading irrelevant perturbation has RG dimen¬ 
sion Ug = —WxY = — 0.785(20)i^ This is related to the 
RG dimension at the decoupled XY FP of the interac¬ 
tion operator between the two complex order parameters. 
Since the RG dimension of the energy operator J d^x (j)^ 
is 1/vxY, the RG dimension of the interspecies interaction 
operator / |(/3ip|(/32p is 

j/„ = 2/i/xy- 3 = -0.0225(4). 


( 11 ) 
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Since Uu < 0, this result implies that the perturbation is 
irrelevant at the decoupled XY FP. However, since = 
—yu is very small, the scaling corrections, that behave as 
_ ^-0.0225^ decay very slowly. 

If the coupling w is negative, only one bosonic compo¬ 
nent is expected to condense. In this case the transition is 
characterized by the symmetry breaking Therefore, 
if the BEC transition is continuous, the critical behav¬ 
ior must belong to another universality class, different 
from the XY one. The corresponding FP has been ex¬ 
tensively studied within the equivalent 0(2)(g)0(2) LGW 
theory4Si^i^ Estimates of the corresponding critical ex¬ 
ponents are: (i) v = 0.57(3) and y = 0.09(1) from 
the resummation of the six-loop expansion within the 
massive zero-momentum schemej^ (ii) v = 0.65(6) and 
r] = 0.09(4) from five-loop calculations within the min¬ 
imal subtraction renormalization schemei^ These theo¬ 
retical results are also supported by experiments, see e.g., 
Ref.i and references thereini^ Therefore, the stable FP 
of the LGW theory ([3]) with attraction domain in the 
region ic < 0 is characterized by the critical exponents 
V Ki 0.6 and rj r; 0.1. Of course, models which are outside 
the attraction domain of the FP are expected to undergo 
a Hrst-order phase transition. 

We finally mention that the critical exponents of the 
unstable 0(4)-symmetric FP along the separatrix w = 0 
are^^ ^ = 0.750(2) and y = 0.0360(3). This FP is 
unstable because the spin-4 perturbation present when 
w ^ 0 has positive RG dimension = 0.125(5) at the 
0(4) FPi ^^i^^ Thus, an 0(4) critical behavior can only be 
observed by performing a proper tuning of the parame¬ 
ters of the model. 

In the following sections we study model (ED in the 
hard-core V —>■ -l-oo limit. Since the intra-species V is 
naively related to the quartic coupling g, we expect g to 
be large in this limit, so that w = g — u > 0. Therefore, 
the BEG transition should be characterized by the simul¬ 
taneous condensation of both components, and controlled 
by the decoupled FP, with a low-temperature phase in 
which both components condense. Scaling corrections, 
due to the on-site density-density interaction between 
the two components, decay very slowly. We also expect 
such corrections to be larger when the inter-species on¬ 
site interaction is attractive, i.e., for U < 0, while they 
should be small in the opposite case U > 0. This is 
also suggested by the fact that the hard-core 2BH model 
becomes equivalent to the one-component model in the 
limit U —>■ - 1 - 00 . Therefore, for U —>■ -l-oo we expect a 
standard XY transition without slowly-decaying 0(^“‘^“) 
scaling corrections. 


C. Multicritical behavior for two unequal bosonic 
gases 

We now discuss a system of two unequal bosonic 
species, such as that described by the more general BH 




FIG. 3: Different phase diagrams for two interacting bosonic 
gases. Thin lines indicate continuous transitions, while the 
thick line represents first-order transitions. RG analyses pre¬ 
dict that only in the tetracritical case (left panel) may one 
have a continuous multicritical behavior at the intersection of 
the transition lines. 

model 

H = ~ ^ ^ ^ ^ i^tx^sv ~ Ms ^ ^ '^SX (12) 

s (xy) S X 

+ ^sxij^sx 1) F G ^ ^ ^lx^2x ■ 

s X X 

In this case we expect a more complex phase diagram, 
showing various phases with transition lines along which 
only one bosonic component condenses, and multicritical 
points (MCPs), where the critical behavior arises from 
the competition of the two distinct U(l) orderings. More 
specifically, a MCP should be observed at the intersection 
of the normal-to-superfluid transition lines where one of 
the components condenses. 

The LGW theory describing the competition of the 
two different U(l) orderings of the model (IT^ is obtained 
by constructing the most general theory of two com¬ 
plex fields (fs{x), with an independent U(l) symmetry for 
each component, without exchange symmetry. It reads 

■Hlgw = j + ^rs\(ps? (13) 

S 

+ ^ + 2w I 

S 

where now we have two quadratic parameters ri and r 2 
and three quartic parameters 51 , 32 , and it. The mul¬ 
ticritical behavior arising from the competition of the 
two distinct U(l) orderings is determined by the RG flow 
when both quadratic parameters ri and r 2 are simulta¬ 
neously tuned to their critical values, keeping the quartic 
parameters gi, g 2 and u fixed. 

The phase diagram of the most general theory, in 
which the associated symmetries are 0 (ni) and 0 ( 112 ), 
has already been investigated within the mean-field 
approximationj^^— Several different phase diagrams 
have been identified, characterized by three or four tran¬ 
sition lines meeting at a MCP, characterized by the pres¬ 
ence or the absence of a mixed phase, in which both fields 
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condense. In Fig. [3] we show the phase diagrams cor¬ 
responding to the case of two coupled U(l)-symnietric 
theories, in the T-S plane where S represents a sec¬ 
ond relevant parameter (for instance, the difference of 
the chemical potentials of the two species) that must be 
tuned to obtain the multicritical behavior. In the LGW 
theory the two behaviors are determined by the sign of 
A = gig2 - u^- If A > 0, four critical lines meet at 
the MCP (tetracritical behavior), as in the left panel of 
Fig. El while, if A < 0, two critical lines and one first- 
order line (bicritical behavior) are present, see the right 
panel of Fig. El 

The sign of A also controls the nature of the behavior 
at the MCP. The FT analysis^i^ of the LGW theory 
0 shows that the system undergoes a first-order tran¬ 
sition at the MCP for A < 0, i.e., in the bicritical case, 
as no stable FP is found in this parameter region. In 
the opposite case (tetracritical phase diagram) a contin¬ 
uous transition is possible at the MCP, controlled by the 
decoupled XY FP. 


III. MEAN-FIELD PHASE DIAGRAM OF THE 
2BH MODEL 

The phase diagram of the 2BH model (P) can be stud¬ 
ied in the mean-field approximation, using 

= [(&L - O + '/'s] [{bsy - (l^s) + <i)s] 
^4>shL+(t>:hsy-\(t>s\\ (14) 

where 4>s = {hsx) are two complex space-independent 
variables, that play the role of order parameters at 
the BEC transitions. Approximation (fldll allows us to 
rewrite Hamiltonian © as a sum of decoupled one-site 
Hamiltonians 

= -2DJY, {cj^shl + M (15) 

S 

-g'^Us + ^"^nsins - I) + Unin2- 

S S 

The spectrum of the theory is invariant under the redef¬ 
inition bs —>■ Ugbs, where Us are two arbitrary phases. 
Therefore, there is no loss of generality if the two param¬ 
eters (ps are assumed to be real. They are determined by 
minimizing the single-site free energy with respect to ps ■ 
In the hard-core limit, since Ug = 0, I, the mean-field 
Hamiltonian (flSl) is defined on a Hilbert space of dimen¬ 
sion 4. In the soft-core case, we may have any occupation 
number, so that the Hilbert space has infinite dimension. 
In practice, we only consider states such that Ug < nmax, 
checking that typical occupation numbers are significanly 
lower than Umax and verifying that the results are stable 
with respect to changes of the cut-off Umax- 

In Fig. P we report \pi\'^ = |(('2p for the hard-core 
model at T = 0 as a function of U and g. The pa¬ 
rameter space is divided into four regions: a central one 



FIG. 4: Zero-temperature U-g phase diagram of the 2BH 
model in the hard-core F —>■ c« limit. 



FIG. 5: Phase diagram of the hard-core 2BH model for 
g = —6, 0, 15 in the U-T plane. In particular, we show the 
normal-to-superfluid transition lines. 



FIG. 6: Phase diagram of the 2BH model for F = 10 and g = 
0 in the U-T plane. The transitions along the line U/V = 1 
separating the two low-temperature phases are of first order. 
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(superfluid domain) in which \(p\'^ ^ 0, and three regions 
in which = 0 that differ by the value of the total 
occupation number n = n\ ni. There is a vacuum re¬ 
gion in which n = 0 and two Mott incompressible phases 
with n = 1 and n = 2. The n = 1 Mott phase presents 
several interesting features related to the isospin degrees 
of freedom per site, which may be described by effective 
low-energy spin Hamiltonians >2^ The superfluid do¬ 
main for p, < 0, separating the n = 0 and n = 2 Mott 
phases, extends in a strip around the line U = 2fi that 
gets narrower as n decreases (the n = 2 Mott phase can¬ 
not include the line U = 2fj, where the on-site interactions 
would cancel for ni = n 2 = !)• This mean-field phase 
diagram appears quite similar to that of the ID Hubbard 
model which can be exactly solved, see, e.g.. Refs. 

(we recall that the ID hard-core 2BH model can be ex¬ 
actly mapped onto the ID fermion Hubbard model). 

The mean-field calculations can be extended to finite 
temperatures by minimizing the free-energy density 

/(</),/?) =-i InZ(</>,/3), Z(<^,/3)=^e-'5^% (16) 

^ i 

with respect to the variational parameters (j)s ■ Ei are the 
energy levels of the single-site Hamiltonian. The finite- 
temperature phase boundaries are obtained by looking 
for the smallest value of /3 = 1 /T for which, at any given 
U and /i, (j)g assume non-zero values. 

The phase diagrams in the hard-core limit are shown 
in Fig. [S] for some values of the chemical potential /r. 
In all cases the low-temperature superfluid phases are 
characterized by the simultaneous condensation of both 
components. In particular, for p. = 0, the case we will 
investigate numerically, there is a low-temperature su¬ 
perfluid phase for U > —6, while the system is always in 
the normal phase for [/ < — 6. 

In Fig. iniwe show the phase diagram for a finite value 
of the intra-species coupling, V — 10, and for p = 0. In 
this case we have two different low-temperature phases: 
when V > U both components condense as in the hard¬ 
core limit, while for V < U only one component con¬ 
denses, breaking the exchange symmetry. These different 
condensed phases are separated by a first-order transi¬ 
tion line at [/ = V. These results are completely con¬ 
sistent with the predictions obtained by analyzing the 
corresponding LGW theory ([3]), see Sec. Ill Al 


IV. NUMERICAL STUDY OF 
TWO-COMPONENT HARD-CORE BOSONS 

We now check some of the theoretical predictions of 
the previous sections. We present a numerical analysis 
of the critical behavior of hard-core 2BH model (H]). As 
discussed in Sec. Ill Al we expect a critical transition in 
the 3D XY universality class with a simultaneous conden¬ 
sation of both components. Correspondingly, we have^ 
V = 0.6717(1), r] = 0.0381(2). However, the asymp¬ 
totic behavior is approached with slowly-decaying scaling 


corrections, that behave as with uju = 0.0225(4). 
These corrections are expected to give rise to significant 
effects when the inter-species on-site interaction is at¬ 
tractive, i.e., for [/ < 0, while they may be negligible in 
the repulsive case. In the following we provide numerical 
evidence for this scenario. 


A. Monte Carlo simulations and observables 


We perform quantum Monte Carlo (QMC) simulations 
of the hard-core 2BH model at zero chemical potential 
/i = 0, on cubic lattices with periodic boundary con¬ 
ditions, for L up to 64. We use the directed operator- 
loop algorithm)^^— which is a particular algorithm using 
the stochastic series expansion (SSE) method.— In the 
simulations we determine the helicity modulus and the 
second-moment correlation length. The helicity modulus 
T is the response of the system to a twist of the boundary 
conditions. It can be obtained from the linear winding 
number Wi along the f**' direction. 


T = 



N+ - N' 


(17) 


where and N~ are the number of non-diagonal oper¬ 
ators which move the particles respectively in the positive 
and negative i^'^ direction. The second-moment correla¬ 
tion length ^ can be conveniently defined from the lattice 
Fourier transform G{p) of the two-point correlation func¬ 
tion G{x — y) = (6^ by), as 


I G(0)^- G{p) 

4sin^(7 r/L) G{p) 


(18) 


where p = {27:jL, 0,0). 

To determine the critical behavior, we perform a finite- 
size scaling (FSS) analysis of the RC invariant quantities 
Ry — TL and — ^/L (we generically denote them as 
R). Close to the transition point {3 = l/T = /3c, they 
behave as 


R{^,L) = /fl(rLl/")+L-“lgfl(TLl/'^)+0(L-“^L-2-l), 

(19) 

where fR{x) is universal apart from a rescaling of its 
argument, r = 1 — (3f Pc, izis the correlation-length expo¬ 
nent, uji (0 < wi < UJ 2 < ...) are the exponents control¬ 
ling the scaling corrections to the asymptotic behavior, 
which are associated with the irrelevant perturbations at 
the stable FP. 

The scaling equation (HU implies that data for dif¬ 
ferent values of L, in particular Li = L and L 2 = 2L, 
cross at a given /3cr(R; Ti, T2), which approaches Pc for 
L —>• 00. More precisely, 

Pcr{R-,L,2L) = Pc + 0{L-^/''-‘^^). ( 20 ) 


Moreover, the value of i? for /3 = Pec approaches the 
universal critical value R* = /fl(0), i.e. 


R{Pcc.L) = R*eY. + b2uL-^^^ + ... ( 21 ) 

n—1 n—1 
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— j^xy+^^xy) 



FIG. 7: (Top) Crossing points Pcr{Rr', L,2L) and 

/?cr(77e; L, 2L) for [7 = 0 and p = 0. They are plotted versus 
I/“1 /‘'xy-“xy _ j^-2.27^ which is the expected behavior of the 
leading scaling corrections. The dashed lines correspond to 
linear fits of the data for the largest available lattices. (Bot¬ 
tom) Estimates of Rr and 7?^ at the crossing points, versus 
_ ^-0.785^ dashed lines correspond to linear fits to 
R* -|-cI/““xY ^ith R* fixed to its 3D XY value = 0.516(1) 
and = 0.5924(4)], which is reported along the ordinate 
axis (crosses). 


B. QMC results for [7 = 0 

To begin with, we consider the hard-core 2BH model 
for U = 0, representing two decoupled and identical 
single-component hard-core BH models. The results will 
then be compared with those obtained for 17 yf 0. 

In this case we have a robust theoretical prediction for 
its critical behavior at the BEG transition: it belongs 
to the 3D XY universality class, described by a stan¬ 
dard U(l)-symmetric theory with one complex order 
parameter i^d^d^ The leading scaling corrections decay 
with exponent wi = ujxy = 0.785(20), and the asymp¬ 
totic critical values of Ry and are Ry = 0.516(1) and 
i?| = 0.5924(4), respectively!^ Numerical evidence of this 
critical behavior has already been reported in Refs. l^ld^ . 

We determine the crossing points^ fici{R] L, 2T) for L 
up to 32. Results are shown in the top panel of Fig. 0 
The behavior for L —^ oo is consistent with the expected 


0 (L“^/‘'xy-<vxy^ scaling corrections. Linear fits of the data 
for L > 20 give /3c = 0.496035(10);^ which is in agree¬ 
ment with, and slightly improves, earlier estimates4^ The 
values of Ry and R^ at the crossing points are reported in 
the lower panel of Fig.[7l They show the expected asymp¬ 
totic behavior R{/3ct,L) = R* + Moreover, the 

extrapolated values are consistent with the best available 
estimates for the 3D XY universality class. For example, 
linear fits of the data for L > 20 give Ry = 0.516(4), 
which is in agreement with the best available estimate^ 
Ry = 0.516(1) of the 3D XY universality class. 

C. QMC results for [7 / 0 

We now present results for the hard-core 2BH model 
for C7 yf 0. In this case we should consider the slowly- 
decaying scaling corrections of order predicted in 

Sec. Ill A[ which are expected to give rise to significant 
systematic deviations, at least for negative U. Since lOu = 
—0.0225, these deviations are hardly detectable numeri¬ 
cally. Indeed, in our range of values of L, 8 < L < 64, 
varies only by 5%, hence it is very difficult to dis¬ 
tinguish it from a constant term. In practice, unless data 
are extremely precise, any FSS analysis is unable to de¬ 
termine the leading scaling function appear¬ 

ing in Eq. m- The extrapolation of the data to L —7 oo 
would identify the asymptotic behavior with that given 

by 

fRirL^/n = fRirL^^n + AgR/rL^I'^), ( 22 ) 

where the slowing-decaying factor is effectively re¬ 
placed with some kind of average A = in the 

considered range of values of L. This observation sug¬ 
gests that the analysis based on the large-L extrapola¬ 
tion of the crossing points should be able to determine 
correctly /3c and v. On the other hand, the extrapolation 
of the data for Ry and at the critical point would 
give /fi(0), which differs from the correct asymptotic es¬ 
timate. In other words, we cannot rely on the values 
of the RG invariant quantities at the critical point to 
identify the universality class. In this discussion we have 
assumed that there is only a single slowly-decaying cor¬ 
rection term, i.e., which is, however, not the case. 

RG also predicts the presence of correction terms pro¬ 
portional to for any integer n, cf. Eq. dUD, which 

makes the analysis of the corrections even more difficult. 
Finally, note that corrections of order , the leading 
ones in the single-species model, are also expected. 

We first consider the attractive hard-core model with 
[/ = —5 and g = 0. Fig. [8] reports the estimates 
of Ry as a function of (3, which clearly show crossing 
points between (3 = 0.714 and /3 = 0.715. The cross¬ 
ing points /3cr(7?T; T, 2L) up to L = 32 are shown in 
the inset. We stress that their determination requires 
no prior knowledge on the nature of the transition, and 
it is therefore completely unbiased. They clearly ap¬ 
pear to converge to a critical value /3c- The precision 







FIG. 8: QMC estimates of Rr for U = —5 and ^ — 0. The 
inset shows the crossing points Pci{Rr', L, 2L) versus . 



FIG. 9: i?T and iZj at the crossing points obtained using 
data for lattices of sizes L and 2L. Here fj, = 0, U = —5 
and U = —3. The straight lines correspond to linear fits to 
7?#(/Icr, L) = with and c as free parameters 

(for U = —5 we obtain RriPa, L) ~ 0.563—0.171/““’''^ for L > 
12 and R^{Pct,L) = 0.618 — 1.401/““’'’' for L > 20). We also 
report the results of the fits to R^{Pci,L) = + c„Z/““" + 

cL-“xy, in which 7?^ is fixed to the 3D XY critical valne and 
Cu and c are free parameters (we obtain Rr{Pci, L) = 0.516 + 
0.052L“““-0.18L““’'’' andi?e(/3cr,7/) = 0.5924+0.0281/“““- 
1.41+““’'’' for U = —5). The crosses on the vertical axis mark 
the best estimates of R* for the 3D universality class. 


of the data does not allow us to distinguish the expected 
(9(2^“1 /i'xy“<v„) approach to the asymptotic value, as pre¬ 
dicted by theory, from the (9(+“1 /'^xy-<^xy^ approach at 
[7 = 0. Nevertheless, we obtain a reasonably precise 
estimate /3c = 0.7144(1), where the error also accommo¬ 
dates the difference of the extrapolations using the two 
ansatzes. Analogous results, although less precise, are 
obtained from the data. 

The values of Rr and at the crossing points are 
shown in Fig. They show an apparently linear be¬ 
havior when plotted against as in the [7 = 0 case. 

However, an extrapolation using the ansatze i?*+cL““’'’', 


which appears consistent with the data, gives critical val¬ 
ues for which are definitely different from 

those of the 3D XY universality class. For example, we 
obtain = 0.563(2) and = 0.625(5) with an accept¬ 
able x^/DOF (DOF is the number of degrees of freedom 
of the fit), which differ significantly from the XY esti¬ 
mates i?Y = 0.516(1) and = 0.5924(4). In view of the 
previous discussion, this discrepancy should be expected, 
because of the presence of the slowly-decaying 0(L“““) 
corrections with lUu = 0.0225(4). For instance, the data 
in Fig. [9] can also be nicely fitted to 

i?(/3„) = i?*+ aL“““+5+““”", (23) 

with R* fixed to its XY value, as shown in Fig. [HI 
Note that fits which include the 0(L“““) corrections 
are hardly distinguishable from the ones that assume the 
leading correction to be in the range of L for which 

the MC data are available. This confirms that disentan¬ 
gling the correction term from the leading constant is ex¬ 
tremely hard, requiring accurate computations for very 
large lattice sizes. It should be stressed that the fit to 
the ansatze (^51) is only an exercise. It is presented to 
make plausible that the transition is in the XY univer¬ 
sality class, even though a naive fit of the data provides 
a different value for R*. Indeed, fit (ESD is not concep¬ 
tually correct unless L is enormously large. RG predicts 
corrections of order for any integer n, that are as 

relevant as the leading one in our range of values of L. 

To further check the above-reported results, we repeat 
the FSS analysis at fixed /3 = 0.7144, varying the on-site 
coupling U which now takes the role that /3 had in the 
previous analysis. An analogous FSS analysis of data up 
to L = 48 gives Uc = —4.9999(3), which is perfectly con¬ 
sistent with the FSS analysis at fixed U. Moreover, the 
values of Rr and R^ at the crossing points in the vari¬ 
able U are hardly distinguishable from those appearing 
in Fig. [5] at the corresponding values of L. 

We have also performed a FSS analysis of data for 
U = —3, up to L = 40, obtaining /3c = 0.5390(2). 
The data of Rr at the crossing points /3cr(7?T; A, 2L) are 
also shown in Fig. [51 As for U = —5, they appear to 
behave linearly with respect to and again they 

extrapolate to a value that is significantly larger than 
Ry Ri 0.516. Such a deviation is smaller than that ob¬ 
tained for U = —5, confirming that the discrepancies 
cannot be interpreted as due to the presence of a new 
universality class. In that case, indeed, one would ob¬ 
tain the same extrapolated value for both values of U. 
Instead, the results are consistent with our RG predic¬ 
tions: the discrepancies increase with |C7|, which is ex¬ 
actly what should be expected if they are related to the 
slowly-decaying corrections due to the interspecies inter¬ 
action. 

We finally mention that, as already anticipated in 
Sec. Ill A[ the scaling corrections induced by the density- 
density on-site interaction turn out to be small when the 
interspecies interaction is repulsive, that is for C7 > 0. 
We have considered two values of [7, G = 1 and 10, and 
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lattices up to L = 40. In both cases, we obtain results 
for Ry and i?| that are in agreement with the XY values. 
Apparently, the slowly-decaying corrections are negligi¬ 
ble, being at most of the size of the statistical errors. 
Note that these corrections vanish for ?7 = 0 (the two 
models are decoupled) and also for U —>■ -boo (the model 
is equivalent to the hard-core BH model for a single bo¬ 
son, hence it has a standard XY transition without the 
scaling corrections). Apparently, they keep on be¬ 
ing small for all intermediate values of U. 

D. Phase diagram of the hard-core 2BH model at 

fj, = 0 

We determine the U dependence of the normal-to- 
superfluid transition line by repeating the FSS analysis 
for other values of U. This is done with less accuracy, us¬ 
ing data up to L = 24. The results are shown in Fig. (TU] 
The phase diagram is quite similar to that obtained in 
the mean-field approximation. In particular, the low- 
temperature superfluid phase disappears for U < — 6 , 
very close to the mean-field result U = — 6. Indeed, as 
shown by the zero-temperature mean-field phase diagram 
shown in Figs. 0] and [Sj C/ = —6 is the location of the 
quantum transition between the superfluid and n = 2 
Mott phase. 

Finally, we discuss the behavior of the particle den¬ 
sities Us = (nsx) at the BEC transition. Their leading 
behavior at the BEC transitions arises from the analyt¬ 
ical background terms, while the universal power laws 
related to the critical behavior are subleading. Indeed, 
standard RG arguments predict the asymptotic behavior 

n««/a(r)+L-^"/,(rLi/") (24) 

when varying the reduced temperature r = 1 — /3//3c 
keeping fixed the model parameters. In Eq. (011), fa is a 
nonuniversal analytic function; is the RG dimension 
of the particle density operator Usx, which is given by 
i/n = ^ — l/v = 1.5112(2) at the decoupled XY FP; /g 
is a universal function apart from a factor and a rescal¬ 
ing of its argument. Therefore, the critical densities at 
Tc are expected to approach a nonuniversal constant in 
the large-L limit. In Fig. [TT] we show the large-L ex¬ 
trapolations of the particle-density data at Tc versus the 
on-site coupling U. The comparison with the mean-field 
computations of Sec. imi shows that the mean-field ap¬ 
proximation of the particle density is quite accurate. 

V. CONCLUSIONS 

We investigate BEC in 3D two-component bosonic sys¬ 
tems. In particular we consider two interacting identical 
bosonic gases, described by the 2BH model (0]), which 
may be interpreted as a lattice two-component bosonic 
system. We study the phase diagram and the critical 
behavior by RG, mean-field, and numerical methods. 


T 



® Fixed U QMC heh 

Fixed T QMC kw 

0 _I_I_I_I_I_I_ 

-6 0 6 12 18 24 
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FIG. 10: Transition line of the hard-core 2BH model for 
= 0, as obtained by the FSS analysis of QMC data. We 
plot Tc = lIPc versus U. Most estimates are obtained in sim¬ 
ulations at fixed U in which T is varied. Close to U = —6, 
where the transition line is almost parallel to the T axis, we 
performed simulations at fixed temperature, determining the 
critical coupling Uc- 



FIG. 11: The single-species particle density ris = {uax) along 
the critical line, as obtained by QMC and mean-field calcula¬ 
tions. 

Our RG analysis is based on a LGW theory with two 
complex scalar fields (associated with the two bosonic 
components), which has the same symmetry as that of 
the bosonic system. In the case of two identical com¬ 
ponents with density-density interactions, the relevant 
global symmetry is 2.2,e ®U(1) 00(1). The mean-field 
analysis predicts two different types of low-temperature 
phases. Depending on the values of the on-site inter¬ 
species U and intra-species V couplings, one may have 
(i) a phase in which the exchange symmetry is conserved 
and both components condense or (ii) a phase in which 
only one component condenses, thus breaking the 22 ,e 
exchange symmetry. 

In case (i), which is generically expected for V > U, 
the transition belongs to the 3D XY universality class. 
More precisely, the critical behavior is controlled by a 
decoupled XY FP, implying an asymptotic decoupling 
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of the critical modes associated with the bosonic compo¬ 
nents. The density-density interaction between these two 
components turns out to be an irrelevant perturbation at 
this FP. It does not affect the asymptotic behavior but 
gives rise to slowly-decaying scaling corrections, that be¬ 
have as where ^ is the diverging length scale at 

transition and a;„ = 3 — 2 /z/xy = 0.0225(4). Of course, 
these slowly-decaying effects are absent at the transition 
of a single bosonic species^^— where the leading scaling 
corrections decay as with Wxy ~ 0.78. The presence 
of slowly decaying corrections makes an accurate check 
of the asymptotic 3D XY critical behavior quite hard, 
essentially because one needs to get very close to the 
critical point to make them negligible. These predictions 
are supported by a FSS analysis of QMC data for the 
2BH model ((T|) in the hard-core limit, i.e., for V —>■ oo. 

The FT RG analysis predicts that the nature of the 
transition should significantly change in the soft-core 
regime, i.e. when V U. In this case only one com¬ 
ponent is expected to condense. The corresponding sym¬ 
metry breaking is therefore different, hence it leads to 
a different universality class in the case of continuous 
transitions. We identify this universality class with that 
of the chiral transition in frustrated two-component spin 
models with noncollinear orderwhich has v si 0.6 
as correlation-length exponent. 

Our RG study also predicts the possibility of a criti¬ 
cal behavior with an extended 0(4) symmetry. However, 
such symmetry enlargement can only be observed by tun¬ 
ing a further parameter beside the temperature. 

It should be stressed that the different critical behav¬ 
iors can only be observed if the system is in the attrac¬ 
tion domain of one of the FPs. If this is not the case, the 
transition would be of first order. 

We also extend our analysis to the more general case 
in which the two bosonic components are not identical. 
The phase diagrams are expected to be more complex, as 
shown in Fig.|3l In particular, they may or may not show 
a low-temperature mixed phase characterized by the con¬ 
densation of both components. According to mean-field 
and RG results, if the mixed phase is present, the phase 
diagram presents a tetracritical point where four tran¬ 
sition lines meet, and the multicritical behavior is con¬ 
trolled by the decoupled XY FP. When the mixed phase 
is absent, thus the low-temperature phases are character¬ 
ized by the BEG of only one component, the transitions 
between the two BEG phases must be first order. In this 
case the competition of the two U(l) orderings does not 
lead to a multicritical behavior, since no stable FPs are 
found in the corresponding parameter region; as a conse¬ 
quence the behavior at the intersection of the transition 
lines is expected to show thermodynamic discontinuities 
analogous to first-order transitions. 

The RG analysis of BEG transitions in mixtures of 
bosonic gases can be straightforwardly extended to the 
two-dimensional (2D) case. In two dimensions bosonic 
systems do not experience BEG. The low-temperature 
phase of single species is characterized by a quasi-long 


range order, where correlations decay algebraically at 
large distances, without the emergence of a nonvanishing 
order parameter. The transitions to this low-temperature 
phase are generally of the Berezinskii-Kosterlitz-Thouless 
(BKT) typej^^— characterized by an exponential in¬ 
crease of the correlation length. The phase diagram of 
mixtures of 2D bosonic systems may show BKT transi¬ 
tions such as those of single bosonic species, and tran¬ 
sitions related to the breaking of the Z2 exchange sym¬ 
metry. A similar situation arises in 2D frustrated two- 
component spin models, see, e.g.. Ref. [t^ and references 
therein. RG scaling arguments analogous to those used 
for the 3D case allow us to infer that the critical be¬ 
havior of identical interacting hard-core bosonic compo¬ 
nents is again controlled by a decoupled BKT FP. Since 
the energy operator is marginal at the BKT transition, 
i.e. Ue = 0, the RG dimension of the density-density 
inter-species coupling is given by ?/„ = 2ye — d = —2 
(d = 2 in this case). Therefore, energy-energy or density- 
density interactions between the bosonic species are ir¬ 
relevant also in two dimensions. Unlike the 3D case, the 
corresponding contributions get rapidly suppressed when 
approaching the critical point, indeed they are 
Therefore, we expect that 2D identical hard-core bosonic 
components, such as the hard-core 2BH model o in 
two dimensions, undergo continuous transitions charac¬ 
terized by decoupled BKT behaviors with multiplicative 
and subleading logarithmic corrections^liZ^ analogous to 
the case of a single bosonic species2&. 

We finally note that cold-atom experiments are usually 
performed in inhomogeneous conditions, due to the pres¬ 
ence of space-dependent trapping forces which effectively 
confine the atomic gas within a limited space regioni^— 
The trapping potential is effectively coupled to the parti¬ 
cle density. Thus, it can be taken into account by adding 
a space-dependent trap term such as 

77trap — 'y ^ U(^)^sa5; (^^) 

X 

to the BH Hamiltonian (El), where Vt is the space- 
dependent potential associated with the external force. 
For example, we may consider Vt{r) = where 

r = I a: I is the distance from the center a; = 0 of the 
trap, which describes a harmonic rotationally-invariant 
trap. The inhomogeneity arising from the trapping po¬ 
tential introduces an additional length scale it into the 
problem, which drastically changes the general features 
of the behavior at a phase transition. For example, the 
correlation functions of the critical modes do not develop 
a diverging length scale in a finite trap. Nevertheless, 
when the trap size it becomes large, we may still observe 
a critical regime around the transition point, with uni¬ 
versal trap-size scaling (TSS) behaviors with respect to 
the trap size it- TSS is controlled by the universality 
class of the phase transition of the homogeneous system. 
It has some analogies with the standard FSS for homoge¬ 
neous systems which we exploited in our numerical study, 
see Sec. El The main difference is that, at the critical 
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point, the correlation length ^ around the center of the 
trap shows a nontrivial power-law dependence on the trap 
size i.e., ^ ^ if where 9 is the universal trap exponent. 
TSS has been numerically checked at the BEC transition 
of a single 3D bosonic gas4^i^ Analogous TSS arguments 
can be applied to the BEC transitions of two-component 
bosonic gases. In particular, the trap exponent in the 
case of the harmonic space-dependence of Vt turns out to 
be21 


9 = 


2u 

l + 2v’ 


(26) 


where v is the correlation-length exponent of the tran¬ 
sition, thus 9 = 0.57327(4) at the decoupled XY fixed 
point controlling the simultaneous condensation of both 
components, and 9 ~ 0.5 when only one component con¬ 
denses. 

Our study is relevant to experiments in which a mix¬ 
ture of two bosonic atomic vapours is cooled to the point 
that at least one of them undergoes a BEC transition and 
the number of particles is separately conserved between 


the two species. Recent years have seen the development 
of many such experiments, with the two bosonic species 
being two hyperfine levels of a single isotope^^— two 
isotopes of the same element^ or heteronuclear mix¬ 
tures of different elementsi^r^. Notably, some mixtures 
were also successfully loaded on optical lattice a^^i^^'^^ . 
The availability of a wide range of atomic species and 
the presence of Eano-Feshbach resonances allow to tune 
the intra-species and inter-species interactions between 
the two components of the mixtures. Additional control 
can be achieved by acting on the depth of an optical lat¬ 
tice. The high degree of tunability of these systems may 
make the direct observation of the transitions we predict 
within the reach of experiments. 
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